model {
for(i in 1:N){
mu_M[i] <- alpha*x[i]
M[i] ~ dnorm(mu_M[i],prec1)

logit(mu_y[i])<- c*x[i]+beta*M[i]
y[i] ~ dbern(mu_y[i])

j1[i]~dunif(0.5,100.5)
j2[i]~dunif(0.5,100.5)
j3[i]<- round(j1[i])
j4[i]<- round(j2[i])

mu_M1[i] <- alpha*(x[i]+deltax)
M1[i] ~ dnorm(mu_M1[i],prec1)

logit(mu_y1[i])<- c*(x[i]+deltax)+beta*M1[i]
te[i]<-(mu_y1[i]-mu_y[i])/deltax

mu_M2[i] <- alpha*(x[j3[i]])
M2[i] ~ dnorm(mu_M2[i],prec1)
logit(mu_y2[i])<- c*x[i]+beta*M2[i]

mu_M3[i] <- alpha*(x[j4[i]])
M3[i] ~ dnorm(mu_M3[i],prec1)
logit(mu_y3[i])<- c*(x[i]+deltax)+beta*M3[i]

de[i]<-(mu_y3[i]-mu_y2[i])/deltax
ie[i]<-te[i]-de[i]
}

alpha ~ dnorm(0.0,1.0E-2)
beta ~ dnorm(0.0,1.0E-2)  
c ~ dnorm(0.0,0.05) 
var1 ~ dgamma(1,0.1)
prec1 <-1/var1
}